# import some utilities
# definition of some plotting tools

%matplotlib inline
from bokeh.layouts import gridplot
from bokeh.layouts import layout
from bokeh.plotting import figure, output_notebook, show
from bokeh.models import ColumnDataSource
from bokeh.plotting import output_file
from bokeh.palettes import Category10
from bokeh.io import export_png
from bokeh.models import Legend
from bokeh.models import Range1d
from bokeh.models import ColorBar
from bokeh.models import LinearAxis
from bokeh.models import DatetimeTickFormatter

# select a palette
#from bokeh.palettes import Dark2_25 as palette
from bokeh.palettes import Spectral11, Category20b
from bokeh.palettes import Spectral6
from bokeh.plotting import output_file
from bokeh.transform import linear_cmap

import itertools  

import datetime

import itertools

import numpy
import os
import re
import sys
from netCDF4 import Dataset

import pandas
import logging
import math

def color_gen():
    yield from itertools.cycle(Category10[10])
color = color_gen()
from bokeh.resources import INLINE
output_notebook(INLINE)
export = False
Loading BokehJS ...
import logging
import sys
sys.path.append(os.path.join(os.environ["HOME"],"Projects","SeaFlect","python"))
from DataPoolAccessModule import AeolusDataClass
from DataPoolAccessModule import QualityControlClass
from SurfaceReflectanceModule import SurfaceReflectanceModelClass


def make_plot(title, hist, edges):
    p = figure(title=title, tools='', background_fill_color="#fafafa")
    p.quad(top=hist, bottom=0, left=edges[:-1], right=edges[1:],
           fill_color="navy", line_color="white", alpha=0.5)
    p.y_range.start = 0
    #p.legend.location = "center_right"
    #p.legend.background_fill_color = "#fefefe"
    p.xaxis.axis_label = 'x'
    p.yaxis.axis_label = 'Pr(x)'
    p.grid.grid_line_color="white"
    return p

verbose_level = 'info'
logger=logging.getLogger()
logging.basicConfig(level=getattr(logging, verbose_level.upper()))
Cm = {}
Li =  {}
wind_speed = numpy.linspace(2, 24, 23)
theory=SurfaceReflectanceModelClass(logger=logger)
theory.Ro_fixed=0.08
Li["low"] = numpy.array([theory.surface_reflectance(model="li",u=key, theta=33) for key in wind_speed]).transpose()    
Cm["low"] = numpy.array([theory.surface_reflectance(model="concept",u=key, theta=33) for key in wind_speed]).transpose()
theory.Ro_fixed=0.016
Li["high"] = numpy.array([theory.surface_reflectance(model="li",u=key, theta=33) for key in wind_speed]).transpose()    
Cm["high"] = numpy.array([theory.surface_reflectance(model="concept",u=key, theta=33) for key in wind_speed]).transpose()
    
AeolusProductObject = AeolusDataClass(debug=False, \
                                          generate=False,
                                          logger=logger)
exp = "fixed_range"
#exp = "flexible_range"

AeolusProductObject.pickle_data_file=os.path.join("SeaFlectProjectData","DataPool","seaflect_data_"+exp + ".pkl")
figure_path = os.path.join("project_figures", "notebook", exp)
if not os.path.exists(figure_path):
    os.makedirs(figure_path)
    
[data_pool, available_periods, available_regions] = AeolusProductObject.get_data()
    
    
AeolusProductObject.data_pool=data_pool
# now that we have a consistent dataset we can make an annual
AeolusProductObject._create_annual_(periods=available_periods)
# and globa
AeolusProductObject._create_global_(regions=available_regions)

#definitions of range bins accoring to c-convention

range_bin={"top":0, "atmosphere":3, "surface":4, "ocean":5, "background":6}

Calibration of return signalΒΆ

The return signal from the range bin is uncalibrated. The Mie useful signal are so-called Acumulated Charge Coupled Device (ACCD) and is a digital number with values depending on the bit-depth of the detector chain. What the calibration correction aims is to convert the ACCD counts to a physical quantity. In this case radiance units for instance $\(mW m^{-1} sr^{-1} cm\)$. The relation we are looking for is like

\[ F^{u} = \alpha \Xi + \beta \]

where \(\alpha\) is the so-called gain, \(\beta\) the offset, \(F^{u}\) and \(\Xi\) the return signal in physical and enginering units (the ACCD-counts) respectively.

The reflected radiation from the range bin is a combination of the return signal from the laser and the reflected solar irradiation

\[ F^{u} = {\cal R}_l F^{d}_l + {\cal R}_s F^{d}_s + \epsilon\]

with \(\epsilon\) an error

It is noted that the first term on the right hand side is a bi-directional term: radiation from a specific direction, is scattered into the \(-180^{o}\) direction, while the second term is a composite of the scattered diffuse and direct radiation from the sun. If we ignore the backscatter from the sun, we are interested in

\[{\cal R}_l = F^{u} / F^{d}_l \]

For the incidence radiation at the top of the range bin \(F^d_l\) we adopt the averaged uv energy from the L1B product stream. This is not correct for the surface range bin, as the laser radiation has to travers through the atmosphere above the surface range bin, so this can only be calulated from a multiple scattering approach. If however we ignore the multiple scattering then we only have to consider the transmission through these atmospheric layers. This is the standard approach. Which for simplicity we ignore here.

The \(F^u\) is estimated from the ACCD counts using the following equation

\[ F^u = K_{mie} * N * c * \Xi \]

where \(K_{mie}\) is a parameter with the same name from the L2A product stream, N the number of pulses within the accumulation period, c a scaling parameter and \(\Xi\) the ACCD-counts. We can combine N and c into one parameter. These terms are needed as we will see that the numerical value of \(K_{mie}\) does not make sense to convert the ACCD counts into physical units. The averaged uv energy has a typical value of 60, \(N = 600 \) and \(K_{mie} = O\left( 10^{12} \right)\). This would imply without an ad-hoc correction, the reflectance from a range bin would be significantly smaller than comparable numbers from independent satellite missions like Sciamachy or Gome where a typical value of the surface reflectance of .02 is found. As we will see the scaling factor aims to bring the Aeolus results within a comparable range of this value and also consistent with model calculations of the surface reflectance.

# energy 
region="E"
period="AnnualCycle"
dataset=data_pool[period][region]
graphic=[]
time = dataset["L1Bray"]["date"]

# Seting the params for the first figure.
p=figure(x_axis_type="datetime", plot_width=1000, plot_height=300, \
         tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
           title="averaged uv energy {} {}".format(period, region))

p.y_range = Range1d(start=55, end=70)
p.extra_y_ranges = {"stdev": Range1d(start=-1, end=1)}

p.circle(time, dataset["L1Bmie"]["avg_uv_energy"][:], color="green", legend_label="Mean")

p.circle(time, dataset["L1Bmie"]["uv_energy_std_dev"][:], color="blue", legend_label="Sdev",\
        y_range_name="stdev")

p.yaxis.axis_label = 'Average UV Energy [-]'
## Adding the second axis to the plot. order is everything  
p.add_layout(LinearAxis(y_range_name="stdev", axis_label='stdev [-]'), 'right')

p.legend.location = "top_left"
graphic.append(p)
plot = gridplot(graphic, ncols=1, width=800, height=300, toolbar_location=None)
show(plot)

if export:
    figure_name = os.path.join(figure_path,"avg_uv_signal_{}_{}.png".format(period, region))
    export_png(plot, filename=figure_name)   
    
graphic=[]   

dataset=data_pool[period][region]
time = dataset["L2A"]["date"]    
i=figure(x_axis_type="datetime",  plot_width=800, plot_height=150, \
         tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
         title="Kmie-Kray {} {}".format(period, region))
i.circle(time, dataset["L2A"]["Kmie"][:]/1.e12, color="green", legend_label="Kmie/1.e12")
i.circle(time, dataset["L2A"]["Kray"][:]/1.e13, color="blue", legend_label="Kray/1.e13")
i.y_range = Range1d(0, 10.)
i.yaxis.axis_label = 'Kmie-Kray [-]'
i.legend.location = "top_left"
graphic.append(i)
plot = gridplot(graphic, ncols=1, width=800, height=300, toolbar_location=None)
show(plot)
if export:
    figure_name = os.path.join(figure_path,"kmie_kray_{}_{}.png".format(period, region))
    export_png(plot, filename=figure_name)   

DiscussionΒΆ

The uv energy has no specific temporal pattern which relates to the previous patterns. this is different for the Kmie which has a clear decreasing trend over time. Why then only the return signal of the upper most range bin has a similar behaviour is not clear. The Kmie and KRay are not understood. And as indicated above the values are too large to be used for a calibration of the ACCD counts. Not shown is that after october 2020 the \(K_{mie}\) values revert to a default value which is off-scale in the above plot.

Proxy calibrated resultsΒΆ

In the next plots the ACCD-counts scaled with the Kmie and an ad-hoc scaling factor are presented. This parameter is referred to as proxy calibrated results, as no radiometric calibration method is available. The adopted scaling values was determined using trial-and-error such that the reflectance of the surface range bin is within an expected range.

region="E"
graphic=[]
period="AnnualCycle"
dataset=data_pool[period][region]
time = dataset["L1Bmie"]["date"]    
t=figure(x_axis_type="datetime", width=800, height=150, \
                  tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
           title="ACCD {} {}".format(period, region))
t.circle(time, dataset["Product"]["calibrated_signal"][:,range_bin["top"]],\
         color="black", legend_label="top")
t.yaxis.axis_label = 'Proxy [-]'
t.y_range = Range1d(0, 20.)
t.legend.location = "top_left"
graphic.append(t)

p=figure(x_axis_type="datetime", width=800, height=150, \
                  tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
           title="ACCD {} {}".format(period, region))
p.circle(time, dataset["Product"]["calibrated_signal"][:,range_bin["atmosphere"]],\
         color="red", legend_label="atmosphere")
p.yaxis.axis_label = 'Proxy  [-]'
p.y_range=Range1d(-10, 20)
p.legend.location = "top_left"
graphic.append(p)

p=figure(x_axis_type="datetime", width=800, height=150, \
                  tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
           title="ACCD {} {}".format(period, region))
p.circle(time, dataset["Product"]["calibrated_signal"][:,range_bin["surface"]], \
         color="green", legend_label="surface")
p.yaxis.axis_label = 'Proxy  [-]'
p.y_range=Range1d(-10, 20)
p.legend.location = "top_left"
graphic.append(p)

p=figure(x_axis_type="datetime", width=800, height=150, \
         tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
         title="ACCD {} {}".format(period, region))
p.circle(time, dataset["Product"]["calibrated_signal"][:,range_bin["ocean"]], \
         color="blue", legend_label="ocean")
p.yaxis.axis_label = 'Proxy  [-]'
p.y_range=Range1d(-20, 20)
p.legend.location = "top_left"
graphic.append(p)

plot = gridplot(graphic, ncols=1, width=800, height=300, toolbar_location="right")
show(plot)
if export:
    figure_name = os.path.join(figure_path,"basic_accd_signal_{}_{}.png".format(period, region))
    export_png(plot, filename=figure_name)   

DiscussionΒΆ

The proxy (as we have no radiometric calibration) return signal for the upper most range bin is decreasing over time. A similar behaviour is not observed in the other signal, in fact there seems to be an increase. At least the dynamical range of the surface range bin (bin 22) tends to be smaller at the end of the time sequence than at the beginning and also the atmospheric range bin (21). These parameters are not corrected for the thickness of the layers, which we know are not constant over time.